#include "PMT.hh"
#include "TFile.h"
#include "TRandom.h"

#include <fstream>
#include <iostream>
#include <string>

Double_t LiheShapeBukin(Double_t* x, Double_t* par)
{
    const double consts = 2 * sqrt(2 * log(2.0));
    double r1 = 0, r2 = 0, r3 = 0, r4 = 0, r5 = 0, hp = 0;
    double x1 = 0, x2 = 0;
    double fit_result = 0;

    double Xp = par[0];
    double sigp = par[1];
    double xi = par[2];
    double rho1 = par[3];
    double rho2 = par[4];

    hp = sigp * consts;
    r3 = log(2.);
    r4 = sqrt(xi * xi + 1);
    r1 = xi / r4;

    if (fabs(xi) > exp(-6.))
    {
        r5 = xi / log(r4 + xi);
    }
    else
        r5 = 1;

    x1 = Xp + (hp / 2) * (r1 - 1);
    x2 = Xp + (hp / 2) * (r1 + 1);

    //--- Left Side
    if (x[0] < x1)
    {
        r2 = rho1 * (x[0] - x1) * (x[0] - x1) / (Xp - x1) / (Xp - x1) - r3 + 4 * r3 * (x[0] - x1) / hp * r5 * r4 / (r4 - xi) / (r4 - xi);
    }

    //--- Center
    else if (x[0] < x2)
    {
        if (fabs(xi) > exp(-6.))
        {
            r2 = log(1 + 4 * xi * r4 * (x[0] - Xp) / hp) / log(1 + 2 * xi * (xi - r4));
            r2 = -r3 * r2 * r2;
        }
        else
        {
            r2 = -4 * r3 * (x[0] - Xp) * (x[0] - Xp) / hp / hp;
        }
    }

    //--- Right Side
    else
    {
        r2 = rho2 * (x[0] - x2) * (x[0] - x2) / (Xp - x2) / (Xp - x2) - r3 - 4 * r3 * (x[0] - x2) / hp * r5 * r4 / (r4 + xi) / (r4 + xi);
    }

    if (fabs(r2) > 100)
    {
        fit_result = 0;
    }
    else
    {
        //---- Normalize the result
        fit_result = exp(r2);
    }

    return fit_result;
}

Double_t LiheShapeNovosibirsks(Double_t* x, Double_t* par)
{
  Double_t peak = par[0];
  Double_t width = par[1];
  Double_t tail = par[2];
  if (TMath::Abs(tail) < 1.e-7) {
    return TMath::Exp( -0.5 * TMath::Power( ( (x[0] - peak) / width ), 2 ));
  }

  Double_t arg = 1.0 - ( x[0] - peak ) * tail / width;

  if (arg < 1.e-7) {
    //Argument of logarithm negative. Real continuation -> function equals zero
    return 0.0;
  }

  Double_t log = TMath::Log(arg);
  static const Double_t xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )

  Double_t width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
  Double_t width_zero2 = width_zero * width_zero;
  Double_t exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );

  return TMath::Exp(exponent) ;

}

PMT::PMT(Int_t iPMT, YAML::Node configer) : m_iPMT(iPMT)
{
    m_VPhoton.clear();
    m_PMTLineShape = nullptr;
    m_SumQ = 0;
    m_TauVal = configer["Tau"].as<double>();
    m_StartTime = configer["StartTime"].as<double>();
    m_StopTime = configer["StopTime"].as<double>();
    YAML::Node configPosition = configer["Position"];
    m_Position = TVector3(configPosition[0].as<double>(),
                          configPosition[1].as<double>(),
                          configPosition[2].as<double>());

    YAML::Node configBaseLine = configer["BaseLineVal"];
    m_BaseLineVal = std::pair<Double_t, Double_t>(configBaseLine["mean"].as<double>(), configBaseLine["sigma"].as<double>());
    YAML::Node configTransit = configer["TransitVal"];
    m_TransitVal = std::pair<Double_t, Double_t>(configTransit["mean"].as<double>(), configTransit["sigma"].as<double>());
    YAML::Node configSmear = configer["SmearVal"];
    m_SmearFunc = new TF1(TString::Format("SmearFunc%d", m_iPMT), LiheShapeBukin, 0.5, 2., 5);
    m_SmearFunc->SetParameters( configSmear["mean"].as<double>(),
                                configSmear["sigma1"].as<double>(), 
                                configSmear["sigma2"].as<double>(), 
                                configSmear["rho1"].as<double>(), 
                                configSmear["rho2"].as<double>());

    YAML::Node configQ = configer["QFunc"];
    m_QFunc = new TF1(TString::Format("QFunc%d", m_iPMT), "TMath::Gaus(x,[0],[1])", 0., 7000.);
    m_QFunc->SetParameters( configQ["mean"].as<double>(),
                            configQ["sigma"].as<double>());

    this->LoadQuantumEff(configer["QuantumEff"].as<std::string>());
}

PMT::~PMT()
{
    this->ClearVPhoton();
    delete m_QFunc;
    delete m_SmearFunc;
    delete m_QuantumEff;
}

void PMT::ClearVPhoton()
{
    for(auto pPtr:m_VPhoton)
    {
        if(pPtr != nullptr)
        {
            delete pPtr;
        }
    }
    m_VPhoton.clear();
    m_SumQ = 0.;    
    if(m_PMTLineShape != nullptr)
    {
        delete m_PMTLineShape;
    }
}

void PMT::ClassifyNewPhoton(aSimPhoton& aPhoton)
{
    if (gRandom->Uniform(100.) < m_QuantumEff->Eval(1242.375/aPhoton.Energy))
    {
        Int_t iPhoton = m_VPhoton.size();
        Double_t Amp = m_QFunc->GetRandom();
        Double_t Time = aPhoton.CatchTime + gRandom->Gaus(m_TransitVal.first, m_TransitVal.second);
        Double_t Smear = m_SmearFunc->GetRandom();
        m_VPhoton.push_back(new PhotonLineShape(m_iPMT, iPhoton, Amp, Time, m_TauVal, Smear, m_StartTime, m_StopTime));
        m_SumQ += Amp;
    }
}

void PMT::LoadQuantumEff(std::string FileName)
{
    std::vector<Double_t> Energy;       Energy.clear();
    std::vector<Double_t> QuantumEff;   QuantumEff.clear();
    std::ifstream InputCSV(FileName.c_str());
    std::string Line;
    while (std::getline(InputCSV, Line))
    {
        std::string Number;
        std::istringstream ReadLine(Line);
        std::pair<Double_t, Double_t> EnergyVSLumi;
        std::getline(ReadLine, Number, '\t');
        Energy.push_back(std::atof(Number.c_str()));
        std::getline(ReadLine, Number, '\t');
        QuantumEff.push_back(std::atof(Number.c_str()));
    }
    m_QuantumEff = new TGraph(Energy.size(), &(Energy[0]), &(QuantumEff[0]));
}

void PMT::Processing()
{
    m_PMTLineShape = nullptr;
    std::vector<Double_t> xVal;     xVal.clear();
    std::vector<Double_t> yVal;     yVal.clear();
    for(double iT = m_StartTime; iT < m_StopTime; ++iT)
    {
        Double_t qVal = gRandom->Gaus(m_BaseLineVal.first, m_BaseLineVal.second);
        for(auto pPtr:m_VPhoton)
        {
            Double_t PhotonV = pPtr->Evaluate(iT);
            qVal = qVal - PhotonV;
        }
        if(qVal < 0)
        {
            qVal = 0.;
        }
        xVal.push_back(iT);
        yVal.push_back(qVal);
    }
    m_PMTLineShape = new TGraph(xVal.size(), &(xVal[0]), &(yVal[0]));
}